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ABSTRACT 


Molecular  dynamics  simulations  which  make  use  of  a  many-body  analytic  potential 
function  have  been  used  to  study  the  nanometer-scale  indentation  of  diamond  and  graphite. 
We  find  that  the  simulation  correctly  reproduces  experimentally  determined  trends  in 
load  versus  penetration  data.  As  a  result,  trends  in  mechanical  properties,  e.g.  Young’s 
modulus,  are  also  reproduced. 


I.  INTRODUCTION 


Microindention  has  long  been  used  to  characterize  bulk  properties  of  materials.  An 
indentation  curve  plots  the  relationship  between  load  and  penetration  depth,  h,  contin¬ 
uously  measured  and  recorded  during  a  hardness  test.1  A  limitation  of  traditional  in¬ 
dentation  techniques  is  the  need  to  image  the  surface  after  indentation  to  determine  the 
penetration  depth  of  the  plastically  deformed  material.  For  microindentation,  SEM  and 
even  TEM  images  are  required.1’2  Furthermore,  shallow  indents  are  often  hard  to  find  and 
difficult  to  image.  As  the  interest  in  thin  films  has  increased,  so  has  the  need  for  reliable 
characterization  methods.  This  and  other  considerations  have  prompted  the  use  of  instru¬ 
ments  such  as  the  Nano-iAdentor3  and  the  atomic  force  microscope  (AFM)2,4,5  to  study 
the  nanomechanical  properties  of  solids  and  thin  films. 

Recently,  the  AFM  has  been  used  to  measure  the  nanomechanical  properties  (e.g. 
elastic  modulus  and  hardness)  of  an  elastomer,  graphite,  gold  films2,  and  diamond 
surfaces,5  and  to  examine  the  surface  forces  of  monolayer  films6  with  depth  and  force 
resolution  superior  to  other  methods.  To  better  understand  the  strengths,  limitations  and 
interpretation  of  nanoindentation  for  characterizing  materials  and  thin  film  properties,  we 
have  been  using  molecular  dynamics  (MD)  to  simulate  indentation  at  the  atomic  scale.  In 
this  work  we  report  initial  studies  of  the  indentation  of  hydrogen  terminated  (111)  and 
(100)  surfaces  of  diamond  and  the  basal  plane  of  graphite  with  a  diamond  tip.  We  find 
that  the  force  on  the  simulated  tip  varies  linearly  with  the  penetration  depth  and  that  the 
slopes  reflect  the  correct  trends  for  the  moduli  of  the  various  surfaces.  We  interpret  the 
differences  in  moduli  for  the  two  diamond  surfaces  in  terms  of  the  relative  displacement  of 
bond  angles  and  bond  lengths  during  indentation. 


II.  APPROACH 


Molecular  dynamics  calculations  were  carried  out  by  integrating  Newton’s  equations 
of  motion  with  a  third  order  variable  time  step  Nordsieck  predictor  corrector.7  The  particle 
forces  are  derived  from  an  empirical  hydrocarbon  potential  that  is  capable  of  modeling 
intramolecular  chemical  bonding  in  both  diamond  and  graphite  lattices,  as  well  as  in  a 
variety  of  small  hydrocarbon  molecules.8  The  potential  used  here  is  potential  II  of  reference 
8  with  additional  terms  that  better  describe  torsional  forces9  and  short  range  repulsive 
terms  which  may  prove  important  under  high  compression.10  This  potential  has  recently 
been  used  to  model  the  compression  of  C$o  between  graphite  planes,9  the  scattering  of  Ceo 
from  diamond  surfaces,11  and  the  compression  and  adhesion  of  diamond  (111)  crystals.10,12 
For  the  graphite  simulations,  interactions  between  graphite  sheets  were  modeled  by  a  pair- 
additive  exponential-6  potential  that  had  been  fit  to  the  experimental  interplanar  distance 
and  compressibility  of  graphite.13 

The  initial  atomic  positions  for  the  tip-substrate  systems  are  shown  in  Fig.  1.  The 
tip  contains  10  layers  of  carbon  atoms  and  two  layers  of  hydrogen  atoms  (325  total  atoms). 
It  was  constructed  by  removing  atoms  from  a  (111)  crystal  until  it  looked  approximately 
like  an  inverted  pyramid  with  a  flattened  apex.  Hydrogen  atoms  were  added  to  the  edge 
carbon  atoms,  i.e.  the  sides  of  the  tip,  to  satisfy  their  valency.  During  the  simulation  the 
top  two  layers  of  the  tip  are  held  rigid,  the  four  layers  closest  to  the  substrate  have  no 
constraints,  and  frictional  forces  are  applied  to  the  middle  six  layers  to  control  the  tem¬ 
perature  of  the  tip.14  The  diamond  (111)  substrate  has  10  carbon  layers  and  two  hydrogen 
layers  (one  on  top  and  one  on  bottom),  each  with  64  atoms  per  layer.  The  initial  surface 
before  indentation  has  the  structure  of  the  truncated  bulk,  in  agreement  with  recent  ex¬ 
perimental  studies.15  The  diamond  (100)  surface  contains  16  carbon  layers,  each  with  64 
atoms  per  layer,  and  one  layer  of  hydrogen  atoms  at  the  top  surface.  At  present  there  is 
no  general  agreement  regarding  the  structure  of  the  (100)  diamond  surface.  Some  experi¬ 
mental  data16  suggest  that  the  bulk  terminated  (1  x  1)  dihydride  (two  hydrogen  atoms  on 
each  carbon  atom)  is  the  most  stable  structure.  By  comparison,  there  are  experimental17 
and  theoretical8,18  data  which  suggest  that  the  dimer  reconstructed  (2x1)  monohydride 
(one  hydrogen  atom  on.  each  carbon  atom)  is  the  most  stable  structure.  Thus,  we  start 
with  the  (2  x  1)  monohydride  structure.  The  graphite  lattice  is  composed  of  five  carbon 
atom  layers  of  128  atoms  each.  The  bottom  two  layers  of  the  diamond  substrates  and  the 
bottom  layer  of  the  graphite  substrate  are  held  rigid.  Frictional  forces  are  applied  to  the 
next  six  layers  of  the  diamond  (111)  substrate,  the  next  seven  layers  of  the  diamond  (100) 
substrate,  and  the  next  four  layers  of  the  graphite  substrate.  The  remaining  top  layers  of 
the  substrates  have  no  constraints.  Periodic  boundary  conditions  are  applied  in  the  x-y 
plane  which  contains  the  (111)  and  (100)  surfaces  of  the  diamond  lattices  and  the  graphite 
sheets  in  order  to  simulate  an  infinite  indentation  surface.  The  temperature  of  the  system 
is  maintained  at  300K  throughout  the  simulation. 

Indentation  is  reformed  in  ffie  followin'!  wav.  The  ririd  lavers  of  the  tip  are  moved 
toward  the  substrate  in  steps  of  0.15  A  .  The  tip-substrate  system  is  then  allowed  to  equi¬ 
librate  for  200-500  time  steps  of  the  order  0.1  femtoseconds.  (Although  these  indentation 
rates  are  orders  of  magnitude  faster  than  typical  experimental  values,  they  are  sufficiently 
slow  to  allow  equilibration  of  the  system.)  Once  the  system  equilibrates,  the  load  on  the 
tip  in  the  z  direction  (parallel  to  the  indentation)  is  calculated  by  summing  the  z  forces 
on  all  the  tip  atoms.  This  number  is  then  averaged  over  the  next  100-250  time  steps.  The 
rigid  layer  of  the  tip  is  then  moved  another  0.15  A  closer  to  the  surface  and  the  above 
process  is  repeated.  Penetration  depth  is  defined  as  the  original  z  position  of  the  substrate 
surface  minus  the  average  z  position  of  the  tip  surface.  In  these  simulations  the  load  is 


A) 


B) 


C) 


Figure  1.  Initial  tip-substrate  configurations:  A)  a  hydrogen  terminated  (lxl)  diamond 
(111)  surface;  B)  a  hydrogen  terminated  (2x1)  diamond  (100)  surface;  and  C) 
the  basal  plane  of  graphite. 

increased  by  decreasing  the  distance  between  the  rigid  layers  of  the  tip  and  the  rigid  layers 
of  the  substrate. 


III.  RESULTS 

Fig.  2  shows  simulated  indentation  (loading)  curves  for  each  of  the  substrates.  The 
curves  are  approximately  linear,  with  slopes  Sm  >  Sioo  >>  S Graphite •  For  the  purposes 
of  this  discussion  only  the  loading  portion  of  the  curves  are  shown.  It  should  be  noted 
that  for  all  cases  examined  here  the  interactions  are  elastic,  i.e.  no  irreparable  damage 
has  occurred  to  the  substrate.  For  comparison  we  show  in  Fig.  3  experimental  indentation 
curves  for  diamond5  and  for  graphite2  measured  with  an  AFM.  (A  tungsten  tip  was  used 
to  indent  the  graphite  and  a  Berkovich  diamond  tip  was  used  to  indent  the  diamond  (100) 
substrate.)  Although  the  scale  of  the  two  sets  of  data  is  much  different,  they  show  similar 
trends  in  that  both  are  approximately  linear,  and  the  curve  for  diamond  has  a  larger  slope 
than  graphite. 

Classical  elasticity  theory  predicts  that  the  slopes  of  the  indentation  curves  should 
be  proportional  to  the  modulus  of  the  material.13  For  cubic  systems,  Young  s  modulus  for 
a  (111)  face,  E(lll),  and  a  (100)  face,  E(100),  is  related  to  the  elastic  constants,  c,j,  of  the 
substrate  via  the  following  relationships20 

E(l  11)  =  6c44(ch  +  2ci2)/(cn  +  2cj2  +  4c44)  (1) 

£(100)  =  cn  +  C]  2  —  2c22/cu.  (2) 

For  hexagonal  substances  such  as  graphite,  Young’s  modulus  is  given  by19 

E  =  [snstn4#  -f  ss3cos40  -(-  (2sj3  +  s44)cos20sm2<?]_1 


(3) 


where  s,j  are  the  compliance  coefficients  and  B  is  the  angle  between  the  long  axis  and  the 
hexagonal  axis.  (See  reference  22  for  the  relationship  between  compliance  coefficients  and 
elastic  constants.)  Substituting  the  elastic  constants  for  the  potential  used  here21  into  the 
Eqs.  1-3  give  values  of  1364  GPa  for  E(lll),  482  GPa  for  E(100)  and  17.3  GPa  ( 0  =  78.7°) 
for  graphite.  Therefore,  in  agreement  with  continuum  theory,  the  differences  in  the  slopes 
of  the  simulated  indentation  curves  yield  the  same  quantitative  trend  as  that  given  using 
the  elastic  constants. 


Figure  2.  Simulated  indentation  curves  for  a  diamond  tip  indenting  graphite  and  diamond. 

Open  diamonds  axe  for  the  hydrogen  terminated  (1  x  1)  diamond  (111)  substrate, 
asterisks  are  for  the  hydrogen  terminated  (2  x  1)  diamond  (100)  substrate,  and 
open  circles  are  for  the  graphite  substrate. 
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Figure  3.  Experimental  load  versus  penetration  depth  for  the  indentation  of  the  basal 
plane  of  graphite  and  the  (100)  surface  of  diamond.  These  data  were  obtained 
with  an  AFM.5 


It  is  encouraging  that  these  simulations  predict  the  proper  trend  in  stiffness  for  the 
diamond  faces.  Because  MD  simulations  allow  the  atomistic  mechanism  of  compression  to 
be  examined,  they  can  yield  insight  into  the  stiffness  difference  between  the  diamond  (111) 
and  (100)  surfaces.  Compare,  for  example,  the  typical  bulk  layer  structure  for  diamond 
(111)  and  diamond  (100)  shown  in  Figure  4.  For  the  (111)  surface,  the  first  layer  atoms  are 
directly  above  the  second  layer  atoms  whereas  for  the  (100)  surface,  the  first  layer  atoms 
are  over  holes  in  the  lattice.  Our  previous  compression  studies  of  diamond  slabs10  showed 


that  it  is  easier  to  bend  a  carbon-carbon  bond  than  to  compress  it.  Therefore,  the  stiffness 
of  the  (111)  surface  is  expected  to  be  greater  than  the  (100)  surface.  Experimentally, 
this  trend  in  Young’s  modulus  is  also  present  (E(lll)=1274  GPa  and  E(1Q0)=1179  GPa), 
however  these  values  are  much  different  than  those  produced  by  our  potential. 


Figure  4.  Side  view  of  the  bulk  lattice  structure  for  A)  for  the  diamond  (111)  lattice  and 
B)  the  diamond  (100)  lattice. 


In  summary,  it  appears  that  molecular  dynamics  simulations  of  indentation  hold 
promise  for  understanding  mechanical  properties  obtained  by  nanoindentation.  The  em¬ 
pirical  potential  used  here  yields  elastic  constants  which  are,  in  general,  too  soft.  Improving 
the  potential  energy  function  so  that  the  elastic  constants  are  better  described  should  allow 
for  a  direct  comparison  of  calculated  mechanical  properties  with  those  obtained  experimen¬ 
tally.  These  simulations  should  also  allow  us  to  compare  simulated  indentation  curves  with 
analytic  models  derived  from  classical  continuum  theory.  Work  in  both  areas  is  currently 
under  way. 
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